#####################################################################################################################

# R CODE FOR THE ARTICLE:
# SHARED POSITIONS ON DIVISIVE BELIEFS EXPLAIN POLICY COLLABORATION: EVIDENCE FROM CLIMATE CHANGE POLICY SUBSYSTEMS IN ELEVEN COUNTRIES"
# CODE BY AASA KARIMO, UNIVERSITY OF HELSINKI & PAUL WAGNER, NORTHUMBRIA UNIVERSITY

#####################################################################################################################

## REQUIRED PACKAGES

if(!require(gdata)) install.packages("gdata",repos = "http://cran.us.r-project.org")
library("gdata")
if(!require(statnet)) install.packages("statnet",repos = "http://cran.us.r-project.org")
library("statnet")
if(!require(texreg)) install.packages("texreg",repos = "http://cran.us.r-project.org")
library("texreg")
if(!require(statnet.common)) install.packages("texreg",repos = "http://cran.us.r-project.org")
library("statnet.common")
if(!require(coda)) install.packages("texreg",repos = "http://cran.us.r-project.org")
library("coda")
if(!require(mice)) install.packages("mice",repos = "http://cran.us.r-project.org")
library("mice")
if(!require(xergm)) install.packages("xergm",repos = "http://cran.us.r-project.org")
library("xergm")
if(!require(ggplot2)) install.packages("ggplot2",repos = "http://cran.us.r-project.org")
library("ggplot2")
if(!require(kableExtra)) install.packages("kableExtra",repos = "http://cran.us.r-project.org")
library("kableExtra")
if(!require(GGally)) install.packages("GGally",repos = "http://cran.us.r-project.org")
library("GGally")
if(!require(scales)) install.packages("scales",repos = "http://cran.us.r-project.org")
library("scales")

#####################################################################################################################

## EXPONENTIAL RANDOM GRAPH MODELS BY COUNTRY

## SET SEED
seed <- 12345
set.seed (seed)

#####################################################################################################################

####################################### AUSTRALIA ###################################################################

# import data files

load("au_n1.RData")
load("au_n3.RData")
load("au_coresim_norm.RData")
load("au_secondarysim_norm.RData")
load("au_divisivesim_norm.RData")
load("au_interact_core_secon.RData")
load("au_interact_core_divisive.RData")
load("au_interact_secon_divisive.RData")

## exponential random graph models

AU_fit1 <- ergm(au_n3 ~ 
                    gwesp(1, fixed = T) +
                    twopath +
                    gwodegree(1, fixed = T) +
                    mutual(same = NULL, diff = FALSE, by = NULL, keep = NULL) +
                    nodeicov("OT1") +
                    nodeicov("n1_indegree") +
                    edgecov(au_coresim_norm) +
                    edgecov(au_n1) +
                    edges, 
                  eval.loglik = TRUE , check.degeneracy = TRUE , 
                  control = control.ergm ( seed = seed , MCMC.samplesize = 10000 , MCMC.interval = 1000))

AU_fit2 <- ergm(au_n3 ~ 
                    gwesp(1, fixed = T) +
                    twopath +
                    gwodegree(1, fixed = T) +
                    mutual(same = NULL, diff = FALSE, by = NULL, keep = NULL) +
                    nodeicov("OT1") +
                    nodeicov("n1_indegree") +
                    edgecov(au_secondarysim_norm) +
                    edgecov(au_n1) +
                    edges, 
                  eval.loglik = TRUE , check.degeneracy = TRUE , 
                  control = control.ergm ( seed = seed , MCMC.samplesize = 10000 , MCMC.interval = 1000))

AU_fit3 <- ergm(au_n3 ~ 
                    gwesp(1, fixed = T) +
                    twopath +
                    gwodegree(1, fixed = T) +
                    mutual(same = NULL, diff = FALSE, by = NULL, keep = NULL) +
                    nodeicov("OT1") +
                    nodeicov("n1_indegree") +
                    edgecov(au_divisivesim_norm) +
                    edgecov(au_n1) +
                    edges, 
                  eval.loglik = TRUE , check.degeneracy = TRUE , 
                  control = control.ergm ( seed = seed , MCMC.samplesize = 10000 , MCMC.interval = 1000))

## exponential random graph models with interactions

AU_fit4 <- ergm(au_n3 ~ 
                    gwesp(1, fixed = T) +
                    twopath +
                    gwodegree(1, fixed = T) +
                    mutual(same = NULL, diff = FALSE, by = NULL, keep = NULL) +
                    nodeicov("OT1") +
                    nodeicov("n1_indegree") +
                    edgecov(au_coresim_norm) +
                    edgecov(au_secondarysim_norm) +
                    edgecov(au_interact_core_secon) +
                    edgecov(au_n1) +
                    edges, 
                  eval.loglik = TRUE , check.degeneracy = TRUE , 
                  control = control.ergm ( seed = seed , MCMC.samplesize = 10000 , MCMC.interval = 1000))

AU_fit5 <- ergm(au_n3 ~ 
                      gwesp(1, fixed = T) +
                      twopath +
                      gwodegree(1, fixed = T) +
                      mutual(same = NULL, diff = FALSE, by = NULL, keep = NULL) +
                      nodeicov("OT1") +
                      nodeicov("n1_indegree") +
                      edgecov(au_coresim_norm) +
                      edgecov(au_divisivesim_norm) +
                      edgecov(au_interact_core_divisive) +
                      edgecov(au_n1) +
                      edges, 
                    eval.loglik = TRUE , check.degeneracy = TRUE , 
                    control = control.ergm ( seed = seed , MCMC.samplesize = 10000 , MCMC.interval = 1000))

AU_fit6 <- ergm(au_n3 ~ 
                      gwesp(1, fixed = T) +
                      twopath +
                      gwodegree(1, fixed = T) +
                      mutual(same = NULL, diff = FALSE, by = NULL, keep = NULL) +
                      nodeicov("OT1") +
                      nodeicov("n1_indegree") +
                      edgecov(au_secondarysim_norm) +
                      edgecov(au_divisivesim_norm) +
                      edgecov(au_interact_secon_divisive) +
                      edgecov(au_n1) +
                      edges, 
                    eval.loglik = TRUE , check.degeneracy = TRUE , 
                    control = control.ergm ( seed = seed , MCMC.samplesize = 10000 , MCMC.interval = 1000))

######################################## BRAZIL ###################################################################

# import data files

load("br_n1.RData")
load("br_n3.RData")
load("br_coresim_norm.RData")
load("br_secondarysim_norm.RData")
load("br_divisivesim_norm.RData")
load("br_interact_core_secon.RData")
load("br_interact_core_divisive.RData")
load("br_interact_secon_divisive.RData")

## exponential random graph models

BR_fit1 <- ergm(br_n3 ~ 
                    gwesp(1, fixed = T) +
                    gwdsp(0.2, fixed = T) +
                    twopath +
                    gwodegree(1, fixed = T) +
                    mutual(same = NULL, diff = FALSE, by = NULL, keep = NULL) +
                    nodeicov("OT1") +
                    nodeicov("n1_indegree") +
                    edgecov(br_coresim_norm) +
                    edgecov(br_n1) +
                    edges, 
                  eval.loglik = TRUE , check.degeneracy = TRUE , 
                  control = control.ergm ( seed = seed , MCMC.samplesize = 10000 , MCMC.interval = 1000))

BR_fit2 <- ergm(br_n3 ~ 
                     gwesp(1, fixed = T) +
                     gwdsp(0.2, fixed = T) +
                     twopath +
                     gwodegree(1, fixed = T) +
                     mutual(same = NULL, diff = FALSE, by = NULL, keep = NULL) +
                     nodeicov("OT1") +
                     nodeicov("n1_indegree") +
                     edgecov(br_secondarysim_norm) +
                     edgecov(br_n1) +
                     edges, 
                   eval.loglik = TRUE , check.degeneracy = TRUE , 
                   control = control.ergm ( seed = seed , MCMC.samplesize = 10000 , MCMC.interval = 1000))

BR_fit3 <- ergm(br_n3 ~ 
                     gwesp(1, fixed = T) +
                     gwdsp(0.2, fixed = T) +
                     twopath +
                     gwodegree(1, fixed = T) +
                     mutual(same = NULL, diff = FALSE, by = NULL, keep = NULL) +
                     nodeicov("OT1") +
                     nodeicov("n1_indegree") +
                     edgecov(br_divisivesim_norm) +
                     edgecov(br_n1) +
                     edges, 
                   eval.loglik = TRUE , check.degeneracy = TRUE , 
                   control = control.ergm ( seed = seed , MCMC.samplesize = 10000 , MCMC.interval = 1000))

## exponential random graph models with interactions

BR_fit4 <- ergm(br_n3 ~ 
                    gwesp(1, fixed = T) +
                    gwdsp(0.2, fixed = T) +
                    twopath +
                    gwodegree(1, fixed = T) +
                    mutual(same = NULL, diff = FALSE, by = NULL, keep = NULL) +
                    nodeicov("OT1") +
                    nodeicov("n1_indegree") +
                    edgecov(br_coresim_norm) +
                    edgecov(br_secondarysim_norm) +
                    edgecov(br_interact_core_secon) +
                    edgecov(br_n1) +
                    edges, 
                  eval.loglik = TRUE , check.degeneracy = TRUE , 
                  control = control.ergm ( seed = seed , MCMC.samplesize = 10000 , MCMC.interval = 1000))

BR_fit5 <- ergm(br_n3 ~ 
                      gwesp(1, fixed = T) +
                      gwdsp(0.2, fixed = T) +
                      twopath +
                      gwodegree(1, fixed = T) +
                      mutual(same = NULL, diff = FALSE, by = NULL, keep = NULL) +
                      nodeicov("OT1") +
                      nodeicov("n1_indegree") +
                      edgecov(br_coresim_norm) +
                      edgecov(br_divisivesim_norm) +
                      edgecov(br_interact_core_divisive) +
                      edgecov(br_n1) +
                      edges, 
                    eval.loglik = TRUE , check.degeneracy = TRUE , 
                    control = control.ergm ( seed = seed , MCMC.samplesize = 10000 , MCMC.interval = 1000))

BR_fit6 <- ergm(br_n3 ~ 
                      gwesp(1, fixed = T) +
                      gwdsp(0.2, fixed = T) +
                      twopath +
                      gwodegree(1, fixed = T) +
                      mutual(same = NULL, diff = FALSE, by = NULL, keep = NULL) +
                      nodeicov("OT1") +
                      nodeicov("n1_indegree") +
                      edgecov(br_secondarysim_norm) +
                      edgecov(br_divisivesim_norm) +
                      edgecov(br_interact_secon_divisive) +
                      edgecov(br_n1) +
                      edges, 
                    eval.loglik = TRUE , check.degeneracy = TRUE , 
                    control = control.ergm ( seed = seed , MCMC.samplesize = 10000 , MCMC.interval = 1000))

####################################### THE CZECH REPUBLIC ################################################################

# import data files

load("cz_n1.RData")
load("cz_n3.RData")
load("cz_coresim_norm.RData")
load("cz_secondarysim_norm.RData")
load("cz_divisivesim_norm.RData")
load("cz_interact_core_secon.RData")
load("cz_interact_core_divisive.RData")
load("cz_interact_secon_divisive.RData")

## exponential random graph models

CZ_fit1 <- ergm(cz_n3 ~ 
                      gwidegree(0.4, fixed = T) +
                      gwesp(1, fixed = T) +
                      gwdsp(0.2, fixed = T) +
                      twopath +
                      gwodegree(1, fixed = T) +
                      mutual(same = NULL, diff = FALSE, by = NULL, keep = NULL) +
                      nodeicov("OT1") +
                      nodeicov("n1_indegree") +
                      edgecov(cz_coresim_norm) +
                      edgecov(cz_n1) +
                      edges, 
                    eval.loglik = TRUE , check.degeneracy = TRUE , 
                    control = control.ergm ( seed = seed , MCMC.samplesize = 10000 , MCMC.interval = 1000))

CZ_fit2 <- ergm(cz_n3 ~ 
                       gwidegree(0.4, fixed = T) +
                       gwesp(1, fixed = T) +
                       gwdsp(0.2, fixed = T) +
                       twopath +
                       gwodegree(1, fixed = T) +
                       mutual(same = NULL, diff = FALSE, by = NULL, keep = NULL) +
                       nodeicov("OT1") +
                       nodeicov("n1_indegree") +
                       edgecov(cz_secondarysim_norm) +
                       edgecov(cz_n1) +
                       edges, 
                     eval.loglik = TRUE , check.degeneracy = TRUE , 
                     control = control.ergm ( seed = seed , MCMC.samplesize = 10000 , MCMC.interval = 1000))

CZ_fit3 <- ergm(cz_n3 ~ 
                       gwidegree(0.4, fixed = T) +
                       gwesp(1, fixed = T) +
                       gwdsp(0.2, fixed = T) +
                       twopath +
                       gwodegree(1, fixed = T) +
                       mutual(same = NULL, diff = FALSE, by = NULL, keep = NULL) +
                       nodeicov("OT1") +
                       nodeicov("n1_indegree") +
                       edgecov(cz_divisivesim_norm) +
                       edgecov(cz_n1) +
                       edges, 
                     eval.loglik = TRUE , check.degeneracy = TRUE , 
                     control = control.ergm ( seed = seed , MCMC.samplesize = 10000 , MCMC.interval = 1000))

## exponential random graph models with interactions

CZ_fit4 <- ergm(cz_n3 ~ 
                      gwidegree(0.4, fixed = T) +
                      gwesp(1, fixed = T) +
                      gwdsp(0.2, fixed = T) +
                      twopath +
                      gwodegree(1, fixed = T) +
                      mutual(same = NULL, diff = FALSE, by = NULL, keep = NULL) +
                      nodeicov("OT1") +
                      nodeicov("n1_indegree") +
                      edgecov(cz_coresim_norm) +
                      edgecov(cz_secondarysim_norm) +
                      edgecov(cz_interact_core_secon) +
                      edgecov(cz_n1) +
                      edges, 
                    eval.loglik = TRUE , check.degeneracy = TRUE , 
                    control = control.ergm ( seed = seed , MCMC.samplesize = 10000 , MCMC.interval = 1000, parallel = 11))

CZ_fit5 <- ergm(cz_n3 ~ 
                       gwidegree(0.4, fixed = T) +
                       gwesp(1, fixed = T) +
                       gwdsp(0.2, fixed = T) +
                       twopath +
                       gwodegree(1, fixed = T) +
                       mutual(same = NULL, diff = FALSE, by = NULL, keep = NULL) +
                       nodeicov("OT1") +
                       nodeicov("n1_indegree") +
                       edgecov(cz_coresim_norm) +
                       edgecov(cz_divisivesim_norm) +
                       edgecov(cz_interact_core_divisive) +
                       edgecov(cz_n1) +
                       edges, 
                     eval.loglik = TRUE , check.degeneracy = TRUE , 
                     control = control.ergm ( seed = seed , MCMC.samplesize = 10000 , MCMC.interval = 1000, parallel = 11))

CZ_fit6 <- ergm(cz_n3 ~ 
                       gwidegree(0.4, fixed = T) +
                       gwesp(1, fixed = T) +
                       gwdsp(0.2, fixed = T) +
                       twopath +
                       gwodegree(1, fixed = T) +
                       mutual(same = NULL, diff = FALSE, by = NULL, keep = NULL) +
                       nodeicov("OT1") +
                       nodeicov("n1_indegree") +
                       edgecov(cz_secondarysim_norm) +
                       edgecov(cz_divisivesim_norm) +
                       edgecov(cz_interact_secon_divisive) +
                       edgecov(cz_n1) +
                       edges, 
                     eval.loglik = TRUE , check.degeneracy = TRUE , 
                     control = control.ergm ( seed = seed , MCMC.samplesize = 10000 , MCMC.interval = 1000, parallel = 11))

############################################# GERMANY #####################################################################

# import data files

load("de_n1.RData")
load("de_n3.RData")
load("de_coresim_norm.RData")
load("de_secondarysim_norm.RData")
load("de_divisivesim_norm.RData")
load("de_interact_core_secon.RData")
load("de_interact_core_divisive.RData")
load("de_interact_secon_divisive.RData")

## exponential random graph models

DE_fit1 <- ergm(de_n3 ~ 
                      gwidegree(0.6, fixed = T) +
                      gwesp(1, fixed = T) +
                      gwdsp(0.1, fixed = T) +
                      twopath +
                      gwodegree(1, fixed = T) +
                      mutual(same = NULL, diff = FALSE, by = NULL, keep = NULL) +
                      nodeicov("OT1") +
                      nodeicov("n1_indegree") +
                      edgecov(de_coresim_norm) +
                      edgecov(de_n1) +
                      edges, 
                    eval.loglik = TRUE , check.degeneracy = TRUE , 
                    control = control.ergm ( seed = seed , MCMC.samplesize = 10000 , MCMC.interval = 1000))

DE_fit2 <- ergm(de_n3 ~ 
                       gwidegree(0.6, fixed = T) +
                       gwesp(1, fixed = T) +
                       gwdsp(0.1, fixed = T) +
                       twopath +
                       gwodegree(1, fixed = T) +
                       mutual(same = NULL, diff = FALSE, by = NULL, keep = NULL) +
                       nodeicov("OT1") +
                       nodeicov("n1_indegree") +
                       edgecov(de_secondarysim_norm) +
                       edgecov(de_n1) +
                       edges, 
                     eval.loglik = TRUE , check.degeneracy = TRUE , 
                     control = control.ergm ( seed = seed , MCMC.samplesize = 10000 , MCMC.interval = 1000))

DE_fit3 <- ergm(de_n3 ~ 
                       gwidegree(0.6, fixed = T) +
                       gwesp(1, fixed = T) +
                       gwdsp(0.1, fixed = T) +
                       twopath +
                       gwodegree(1, fixed = T) +
                       mutual(same = NULL, diff = FALSE, by = NULL, keep = NULL) +
                       nodeicov("OT1") +
                       nodeicov("n1_indegree") +
                       edgecov(de_divisivesim_norm) +
                       edgecov(de_n1) +
                       edges, 
                     eval.loglik = TRUE , check.degeneracy = TRUE , 
                     control = control.ergm ( seed = seed , MCMC.samplesize = 10000 , MCMC.interval = 1000))

## exponential random graph models with interactions

DE_fit4 <- ergm(de_n3 ~ 
                        gwidegree(0.6, fixed = T) +
                        gwesp(1, fixed = T) +
                        gwdsp(0.1, fixed = T) +
                        twopath +
                        gwodegree(1, fixed = T) +
                        mutual(same = NULL, diff = FALSE, by = NULL, keep = NULL) +
                        nodeicov("OT1") +
                        nodeicov("n1_indegree") +
                        edgecov(de_coresim_norm) +
                        edgecov(de_secondarysim_norm) +
                        edgecov(de_interact_core_secon) +
                        edgecov(de_n1) +
                        edges, 
                      eval.loglik = TRUE , check.degeneracy = TRUE , 
                      control = control.ergm ( seed = seed , MCMC.samplesize = 10000 , MCMC.interval = 1000, parallel = 11))

DE_fit5 <- ergm(de_n3 ~ 
                        gwidegree(0.6, fixed = T) +
                        gwesp(1, fixed = T) +
                        gwdsp(0.1, fixed = T) +
                        twopath +
                        gwodegree(1, fixed = T) +
                        mutual(same = NULL, diff = FALSE, by = NULL, keep = NULL) +
                        nodeicov("OT1") +
                        nodeicov("n1_indegree") +
                        edgecov(de_coresim_norm) +
                        edgecov(de_divisivesim_norm) +
                        edgecov(de_interact_core_divisive) +
                        edgecov(de_n1) +
                        edges, 
                      eval.loglik = TRUE , check.degeneracy = TRUE , 
                      control = control.ergm ( seed = seed , MCMC.samplesize = 10000 , MCMC.interval = 1000, parallel = 11))

DE_fit6 <- ergm(de_n3 ~ 
                        gwidegree(0.6, fixed = T) +
                        gwesp(1, fixed = T) +
                        gwdsp(0.1, fixed = T) +
                        twopath +
                        gwodegree(1, fixed = T) +
                        mutual(same = NULL, diff = FALSE, by = NULL, keep = NULL) +
                        nodeicov("OT1") +
                        nodeicov("n1_indegree") +
                        edgecov(de_secondarysim_norm) +
                        edgecov(de_divisivesim_norm) +
                        edgecov(de_interact_secon_divisive) +
                        edgecov(de_n1) +
                        edges, 
                      eval.loglik = TRUE , check.degeneracy = TRUE , 
                      control = control.ergm ( seed = seed , MCMC.samplesize = 10000 , MCMC.interval = 1000, parallel = 11))

############################################# FINLAND ######################################################################

# import data files

load("fi_n1.RData")
load("fi_n3.RData")
load("fi_coresim_norm.RData")
load("fi_secondarysim_norm.RData")
load("fi_divisivesim_norm.RData")
load("fi_interact_core_secon.RData")
load("fi_interact_core_divisive.RData")
load("fi_interact_secon_divisive.RData")

## exponential random graph models

FI_fit1 <- ergm(fi_n3 ~ 
                    gwesp(1, fixed = T) +
                    gwdsp(0.1, fixed = T) +
                    twopath +
                    gwodegree(3, fixed = T) +
                    mutual(same = NULL, diff = FALSE, by = NULL, keep = NULL) +
                    nodeicov("OT1") +
                    nodeicov("n1_indegree") +
                    edgecov(fi_coresim_norm) +
                    edgecov(fi_n1) +
                    edges, 
                  eval.loglik = TRUE , check.degeneracy = TRUE , 
                  control = control.ergm ( seed = seed , MCMC.samplesize = 10000 , MCMC.interval = 1000))

FI_fit2 <- ergm(fi_n3 ~ 
                     gwesp(1, fixed = T) +
                     gwdsp(0.1, fixed = T) +
                     twopath +
                     gwodegree(3, fixed = T) +
                     mutual(same = NULL, diff = FALSE, by = NULL, keep = NULL) +
                     nodeicov("OT1") +
                     nodeicov("n1_indegree") +
                     edgecov(fi_secondarysim_norm) +
                     edgecov(fi_n1) +
                     edges, 
                   eval.loglik = TRUE , check.degeneracy = TRUE , 
                   control = control.ergm ( seed = seed , MCMC.samplesize = 10000 , MCMC.interval = 1000))

FI_fit3 <- ergm(fi_n3 ~ 
                     gwesp(1, fixed = T) +
                     gwdsp(0.1, fixed = T) +
                     twopath +
                     gwodegree(3, fixed = T) +
                     mutual(same = NULL, diff = FALSE, by = NULL, keep = NULL) +
                     nodeicov("OT1") +
                     nodeicov("n1_indegree") +
                     edgecov(fi_divisivesim_norm) +
                     edgecov(fi_n1) +
                     edges, 
                   eval.loglik = TRUE , check.degeneracy = TRUE , 
                   control = control.ergm ( seed = seed , MCMC.samplesize = 10000 , MCMC.interval = 1000))

## exponential random graph models with interactions

FI_fit4 <- ergm(fi_n3 ~ 
                      gwesp(1, fixed = T) +
                      gwdsp(0.1, fixed = T) +
                      twopath +
                      gwodegree(3, fixed = T) +
                      mutual(same = NULL, diff = FALSE, by = NULL, keep = NULL) +
                      nodeicov("OT1") +
                      nodeicov("n1_indegree") +
                      edgecov(fi_coresim_norm) +
                      edgecov(fi_secondarysim_norm) +
                      edgecov(fi_interact_core_secon) +
                      edgecov(fi_n1) +
                      edges, 
                    eval.loglik = TRUE , check.degeneracy = TRUE , 
                    control = control.ergm ( seed = seed , MCMC.samplesize = 10000 , MCMC.interval = 1000, parallel = 11))

FI_fit5 <- ergm(fi_n3 ~ 
                      gwesp(1, fixed = T) +
                      gwdsp(0.1, fixed = T) +
                      twopath +
                      gwodegree(3, fixed = T) +
                      mutual(same = NULL, diff = FALSE, by = NULL, keep = NULL) +
                      nodeicov("OT1") +
                      nodeicov("n1_indegree") +
                      edgecov(fi_coresim_norm) +
                      edgecov(fi_divisivesim_norm) +
                      edgecov(fi_interact_core_divisive) +
                      edgecov(fi_n1) +
                      edges, 
                    eval.loglik = TRUE , check.degeneracy = TRUE , 
                    control = control.ergm ( seed = seed , MCMC.samplesize = 10000 , MCMC.interval = 1000, parallel = 11))

FI_fit6 <- ergm(fi_n3 ~ 
                      gwesp(1, fixed = T) +
                      gwdsp(0.1, fixed = T) +
                      twopath +
                      gwodegree(3, fixed = T) +
                      mutual(same = NULL, diff = FALSE, by = NULL, keep = NULL) +
                      nodeicov("OT1") +
                      nodeicov("n1_indegree") +
                      edgecov(fi_secondarysim_norm) +
                      edgecov(fi_divisivesim_norm) +
                      edgecov(fi_interact_secon_divisive) +
                      edgecov(fi_n1) +
                      edges, 
                    eval.loglik = TRUE , check.degeneracy = TRUE , 
                    control = control.ergm ( seed = seed , MCMC.samplesize = 10000 , MCMC.interval = 1000, parallel = 11))

################################################ IRELAND #################################################################

# import data files

load("ie_n1.RData")
load("ie_n3.RData")
load("ie_coresim_norm.RData")
load("ie_secondarysim_norm.RData")
load("ie_divisivesim_norm.RData")
load("ie_interact_core_secon.RData")
load("ie_interact_core_divisive.RData")
load("ie_interact_secon_divisive.RData")

## exponential random graph models

IE_fit1 <- ergm(ie_n3 ~ 
                      gwidegree(0.4, fixed = T) +
                      gwesp(1, fixed = T) +
                      twopath +
                      gwodegree(0.9, fixed = T) +
                      mutual(same = NULL, diff = FALSE, by = NULL, keep = NULL) +
                      nodeicov("OT1") +
                      nodeicov("n1_indegree") +
                      edgecov(ie_coresim_norm) +
                      edgecov(ie_n1) +
                      edges, 
                    eval.loglik = TRUE , check.degeneracy = TRUE , 
                    control = control.ergm ( seed = seed , MCMC.samplesize = 10000 , MCMC.interval = 1000))

IE_fit2 <- ergm(ie_n3 ~ 
                       gwidegree(0.4, fixed = T) +
                       gwesp(1, fixed = T) +
                       twopath +
                       gwodegree(0.9, fixed = T) +
                       mutual(same = NULL, diff = FALSE, by = NULL, keep = NULL) +
                       nodeicov("OT1") +
                       nodeicov("n1_indegree") +
                       edgecov(ie_secondarysim_norm) +
                       edgecov(ie_n1) +
                       edges, 
                     eval.loglik = TRUE , check.degeneracy = TRUE , 
                     control = control.ergm ( seed = seed , MCMC.samplesize = 10000 , MCMC.interval = 1000))

IE_fit3 <- ergm(ie_n3 ~ 
                       gwidegree(0.4, fixed = T) +
                       gwesp(1, fixed = T) +
                       twopath +
                       gwodegree(0.9, fixed = T) +
                       mutual(same = NULL, diff = FALSE, by = NULL, keep = NULL) +
                       nodeicov("OT1") +
                       nodeicov("n1_indegree") +
                       edgecov(ie_divisivesim_norm) +
                       edgecov(ie_n1) +
                       edges, 
                     eval.loglik = TRUE , check.degeneracy = TRUE , 
                     control = control.ergm ( seed = seed , MCMC.samplesize = 10000 , MCMC.interval = 1000))

## exponential random graph models with interactions

IE_fit4 <- ergm(ie_n3 ~ 
                        gwidegree(0.4, fixed = T) +
                        gwesp(1, fixed = T) +
                        twopath +
                        gwodegree(0.9, fixed = T) +
                        mutual(same = NULL, diff = FALSE, by = NULL, keep = NULL) +
                        nodeicov("OT1") +
                        nodeicov("n1_indegree") +
                        edgecov(ie_coresim_norm) +
                        edgecov(ie_secondarysim_norm) +
                        edgecov(ie_interact_core_secon) +
                        edgecov(ie_n1) +
                        edges, 
                      eval.loglik = TRUE , check.degeneracy = TRUE , 
                      control = control.ergm ( seed = seed , MCMC.samplesize = 10000 , MCMC.interval = 1000, parallel = 11))

IE_fit5 <- ergm(ie_n3 ~ 
                        gwidegree(0.4, fixed = T) +
                        gwesp(1, fixed = T) +
                        twopath +
                        gwodegree(0.9, fixed = T) +
                        mutual(same = NULL, diff = FALSE, by = NULL, keep = NULL) +
                        nodeicov("OT1") +
                        nodeicov("n1_indegree") +
                        edgecov(ie_coresim_norm) +
                        edgecov(ie_divisivesim_norm) +
                        edgecov(ie_interact_core_divisive) +
                        edgecov(ie_n1) +
                        edges, 
                      eval.loglik = TRUE , check.degeneracy = TRUE , 
                      control = control.ergm ( seed = seed , MCMC.samplesize = 10000 , MCMC.interval = 1000, parallel = 11))

IE_fit6 <- ergm(ie_n3 ~ 
                        gwidegree(0.4, fixed = T) +
                        gwesp(1, fixed = T) +
                        twopath +
                        gwodegree(0.9, fixed = T) +
                        mutual(same = NULL, diff = FALSE, by = NULL, keep = NULL) +
                        nodeicov("OT1") +
                        nodeicov("n1_indegree") +
                        edgecov(ie_secondarysim_norm) +
                        edgecov(ie_divisivesim_norm) +
                        edgecov(ie_interact_secon_divisive) +
                        edgecov(ie_n1) +
                        edges, 
                      eval.loglik = TRUE , check.degeneracy = TRUE , 
                      control = control.ergm ( seed = seed , MCMC.samplesize = 10000 , MCMC.interval = 1000, parallel = 11))

############################################### JAPAN ######################################################################

# import data files

load("jp_n1.RData")
load("jp_n3.RData")
load("jp_coresim_norm.RData")
load("jp_secondarysim_norm.RData")
load("jp_divisivesim_norm.RData")
load("jp_interact_core_secon.RData")
load("jp_interact_core_divisive.RData")
load("jp_interact_secon_divisive.RData")

## exponential random graph models

JP_fit1 <- ergm(jp_n3 ~ 
                    gwesp(1, fixed = T) +
                    twopath +
                    gwodegree(0.9, fixed = T) +
                    mutual(same = NULL, diff = FALSE, by = NULL, keep = NULL) +
                    nodeicov("OT1") +
                    nodeicov("n1_indegree") +
                    edgecov(jp_coresim_norm) +
                    edgecov(jp_n1) +
                    edges, 
                  eval.loglik = TRUE , check.degeneracy = TRUE , 
                  control = control.ergm ( seed = seed , MCMC.samplesize = 10000 , MCMC.interval = 1000))

JP_fit2 <- ergm(jp_n3 ~ 
                     gwesp(1, fixed = T) +
                     twopath +
                     gwodegree(0.9, fixed = T) +
                     mutual(same = NULL, diff = FALSE, by = NULL, keep = NULL) +
                     nodeicov("OT1") +
                     nodeicov("n1_indegree") +
                     edgecov(jp_secondarysim_norm) +
                     edgecov(jp_n1) +
                     edges, 
                   eval.loglik = TRUE , check.degeneracy = TRUE , 
                   control = control.ergm ( seed = seed , MCMC.samplesize = 10000 , MCMC.interval = 1000))

JP_fit3 <- ergm(jp_n3 ~ 
                     gwesp(1, fixed = T) +
                     twopath +
                     gwodegree(0.9, fixed = T) +
                     mutual(same = NULL, diff = FALSE, by = NULL, keep = NULL) +
                     nodeicov("OT1") +
                     nodeicov("n1_indegree") +
                     edgecov(jp_divisivesim_norm) +
                     edgecov(jp_n1) +
                     edges, 
                   eval.loglik = TRUE , check.degeneracy = TRUE , 
                   control = control.ergm ( seed = seed , MCMC.samplesize = 10000 , MCMC.interval = 1000))

## exponential random graph models with interactions

JP_fit4 <- ergm(jp_n3 ~ 
                      gwesp(1, fixed = T) +
                      twopath +
                      gwodegree(0.9, fixed = T) +
                      mutual(same = NULL, diff = FALSE, by = NULL, keep = NULL) +
                      nodeicov("OT1") +
                      nodeicov("n1_indegree") +
                      edgecov(jp_coresim_norm) +
                      edgecov(jp_secondarysim_norm) +
                      edgecov(jp_interact_core_secon) +
                      edgecov(jp_n1) +
                      edges, 
                    eval.loglik = TRUE , check.degeneracy = TRUE , 
                    control = control.ergm ( seed = seed , MCMC.samplesize = 10000 , MCMC.interval = 1000, parallel = 11))

JP_fit5 <- ergm(jp_n3 ~ 
                      gwesp(1, fixed = T) +
                      twopath +
                      gwodegree(0.9, fixed = T) +
                      mutual(same = NULL, diff = FALSE, by = NULL, keep = NULL) +
                      nodeicov("OT1") +
                      nodeicov("n1_indegree") +
                      edgecov(jp_coresim_norm) +
                      edgecov(jp_divisivesim_norm) +
                      edgecov(jp_interact_core_divisive) +
                      edgecov(jp_n1) +
                      edges, 
                    eval.loglik = TRUE , check.degeneracy = TRUE , 
                    control = control.ergm ( seed = seed , MCMC.samplesize = 10000 , MCMC.interval = 1000, parallel = 11))

JP_fit6 <- ergm(jp_n3 ~ 
                      gwesp(1, fixed = T) +
                      twopath +
                      gwodegree(0.9, fixed = T) +
                      mutual(same = NULL, diff = FALSE, by = NULL, keep = NULL) +
                      nodeicov("OT1") +
                      nodeicov("n1_indegree") +
                      edgecov(jp_secondarysim_norm) +
                      edgecov(jp_divisivesim_norm) +
                      edgecov(jp_interact_secon_divisive) +
                      edgecov(jp_n1) +
                      edges, 
                    eval.loglik = TRUE , check.degeneracy = TRUE , 
                    control = control.ergm ( seed = seed , MCMC.samplesize = 10000 , MCMC.interval = 1000, parallel = 11))

############################################# KOREA ######################################################################

# import data files

load("kr_n1.RData")
load("kr_n3.RData")
load("kr_coresim_norm.RData")
load("kr_secondarysim_norm.RData")
load("kr_divisivesim_norm.RData")
load("kr_interact_core_secon.RData")
load("kr_interact_core_divisive.RData")
load("kr_interact_secon_divisive.RData")

## exponential random graph models

KR_fit1 <- ergm(kr_n3 ~ 
                      gwidegree(0.9, fixed =T) +
                      gwesp(1, fixed = T) +
                      gwdsp(0.1, fixed = T) +
                      twopath +
                      gwodegree(1, fixed = T) +
                      mutual(same = NULL, diff = FALSE, by = NULL, keep = NULL) +
                      nodeicov("OT1") +
                      nodeicov("n1_indegree") +
                      edgecov(kr_coresim_norm) +
                      edgecov(kr_n1) +
                      edges, 
                    eval.loglik = TRUE , check.degeneracy = TRUE , 
                    control = control.ergm ( seed = seed , MCMC.samplesize = 10000 , MCMC.interval = 1000))

KR_fit2 <- ergm(kr_n3 ~ 
                       gwidegree(0.9, fixed =T) +
                       gwesp(1, fixed = T) +
                       gwdsp(0.1, fixed = T) +
                       twopath +
                       gwodegree(1, fixed = T) +
                       mutual(same = NULL, diff = FALSE, by = NULL, keep = NULL) +
                       nodeicov("OT1") +
                       nodeicov("n1_indegree") +
                       edgecov(kr_secondarysim_norm) +
                       edgecov(kr_n1) +
                       edges, 
                     eval.loglik = TRUE , check.degeneracy = TRUE , 
                     control = control.ergm ( seed = seed , MCMC.samplesize = 10000 , MCMC.interval = 1000))

KR_fit3 <- ergm(kr_n3 ~ 
                       gwidegree(0.9, fixed =T) +
                       gwesp(1, fixed = T) +
                       gwdsp(0.1, fixed = T) +
                       twopath +
                       gwodegree(1, fixed = T) +
                       mutual(same = NULL, diff = FALSE, by = NULL, keep = NULL) +
                       nodeicov("OT1") +
                       nodeicov("n1_indegree") +
                       edgecov(kr_divisivesim_norm) +
                       edgecov(kr_n1) +
                       edges, 
                     eval.loglik = TRUE , check.degeneracy = TRUE , 
                     control = control.ergm ( seed = seed , MCMC.samplesize = 10000 , MCMC.interval = 1000))

## exponential random graph models with interactions

KR_fit4 <- ergm(kr_n3 ~ 
                        gwidegree(0.9, fixed =T) +
                        gwesp(1, fixed = T) +
                        gwdsp(0.1, fixed = T) +
                        twopath +
                        gwodegree(1, fixed = T) +
                        mutual(same = NULL, diff = FALSE, by = NULL, keep = NULL) +
                        nodeicov("OT1") +
                        nodeicov("n1_indegree") +
                        edgecov(kr_coresim_norm) +
                        edgecov(kr_secondarysim_norm) +
                        edgecov(kr_interact_core_secon) +
                        edgecov(kr_n1) +
                        edges, 
                      eval.loglik = TRUE , check.degeneracy = TRUE , 
                      control = control.ergm ( seed = seed , MCMC.samplesize = 10000 , MCMC.interval = 1000, parallel = 11))

KR_fit5 <- ergm(kr_n3 ~ 
                        gwidegree(0.9, fixed =T) +
                        gwesp(1, fixed = T) +
                        gwdsp(0.1, fixed = T) +
                        twopath +
                        gwodegree(1, fixed = T) +
                        mutual(same = NULL, diff = FALSE, by = NULL, keep = NULL) +
                        nodeicov("OT1") +
                        nodeicov("n1_indegree") +
                        edgecov(kr_coresim_norm) +
                        edgecov(kr_divisivesim_norm) +
                        edgecov(kr_interact_core_divisive) +
                        edgecov(kr_n1) +
                        edges, 
                      eval.loglik = TRUE , check.degeneracy = TRUE , 
                      control = control.ergm ( seed = seed , MCMC.samplesize = 10000 , MCMC.interval = 1000, parallel = 11))

KR_fit6 <- ergm(kr_n3 ~ 
                        gwidegree(0.9, fixed =T) +
                        gwesp(1, fixed = T) +
                        gwdsp(0.1, fixed = T) +
                        twopath +
                        gwodegree(1, fixed = T) +
                        mutual(same = NULL, diff = FALSE, by = NULL, keep = NULL) +
                        nodeicov("OT1") +
                        nodeicov("n1_indegree") +
                        edgecov(kr_secondarysim_norm) +
                        edgecov(kr_divisivesim_norm) +
                        edgecov(kr_interact_secon_divisive) +
                        edgecov(kr_n1) +
                        edges, 
                      eval.loglik = TRUE , check.degeneracy = TRUE , 
                      control = control.ergm ( seed = seed , MCMC.samplesize = 10000 , MCMC.interval = 1000, parallel = 11))

############################################# PORTUGAL #####################################################################
# import data files

load("pt_n1.RData")
load("pt_n3.RData")
load("pt_coresim_norm.RData")
load("pt_secondarysim_norm.RData")
load("pt_divisivesim_norm.RData")
load("pt_interact_core_secon.RData")
load("pt_interact_core_divisive.RData")
load("pt_interact_secon_divisive.RData")

## exponential random graph models

PT_fit1 <- ergm(pt_n3 ~ 
                    gwesp(1, fixed = T) +
                    twopath +
                    gwodegree(0.8, fixed = T) +
                    mutual(same = NULL, diff = FALSE, by = NULL, keep = NULL) +
                    nodeicov("OT1") +
                    nodeicov("n1_indegree") +
                    edgecov(pt_coresim_norm) +
                    edgecov(pt_n1) +
                    edges, 
                  eval.loglik = TRUE , check.degeneracy = TRUE , 
                  control = control.ergm ( seed = seed , MCMC.samplesize = 10000 , MCMC.interval = 1000))

PT_fit2 <- ergm(pt_n3 ~ 
                     gwesp(1, fixed = T) +
                     twopath +
                     gwodegree(0.8, fixed = T) +
                     mutual(same = NULL, diff = FALSE, by = NULL, keep = NULL) +
                     nodeicov("OT1") +
                     nodeicov("n1_indegree") +
                     edgecov(pt_secondarysim_norm) +
                     edgecov(pt_n1) +
                     edges, 
                   eval.loglik = TRUE , check.degeneracy = TRUE , 
                   control = control.ergm ( seed = seed , MCMC.samplesize = 10000 , MCMC.interval = 1000))

PT_fit3 <- ergm(pt_n3 ~ 
                     gwesp(1, fixed = T) +
                     twopath +
                     gwodegree(0.8, fixed = T) +
                     mutual(same = NULL, diff = FALSE, by = NULL, keep = NULL) +
                     nodeicov("OT1") +
                     nodeicov("n1_indegree") +
                     edgecov(pt_divisivesim_norm) +
                     edgecov(pt_n1) +
                     edges, 
                   eval.loglik = TRUE , check.degeneracy = TRUE , 
                   control = control.ergm ( seed = seed , MCMC.samplesize = 10000 , MCMC.interval = 1000))

## exponential random graph models with interactions

PT_fit4 <- ergm(pt_n3 ~ 
                      gwesp(1, fixed = T) +
                      twopath +
                      gwodegree(0.8, fixed = T) +
                      mutual(same = NULL, diff = FALSE, by = NULL, keep = NULL) +
                      nodeicov("OT1") +
                      nodeicov("n1_indegree") +
                      edgecov(pt_coresim_norm) +
                      edgecov(pt_secondarysim_norm) +
                      edgecov(pt_interact_core_secon) +
                      edgecov(pt_n1) +
                      edges, 
                    eval.loglik = TRUE , check.degeneracy = TRUE , 
                    control = control.ergm ( seed = seed , MCMC.samplesize = 10000 , MCMC.interval = 1000, parallel = 11))

PT_fit5 <- ergm(pt_n3 ~ 
                      gwesp(1, fixed = T) +
                      twopath +
                      gwodegree(0.8, fixed = T) +
                      mutual(same = NULL, diff = FALSE, by = NULL, keep = NULL) +
                      nodeicov("OT1") +
                      nodeicov("n1_indegree") +
                      edgecov(pt_coresim_norm) +
                      edgecov(pt_divisivesim_norm) +
                      edgecov(pt_interact_core_divisive) +
                      edgecov(pt_n1) +
                      edges, 
                    eval.loglik = TRUE , check.degeneracy = TRUE , 
                    control = control.ergm ( seed = seed , MCMC.samplesize = 10000 , MCMC.interval = 1000, parallel = 11))

PT_fit6 <- ergm(pt_n3 ~ 
                      gwesp(1, fixed = T) +
                      twopath +
                      gwodegree(0.8, fixed = T) +
                      mutual(same = NULL, diff = FALSE, by = NULL, keep = NULL) +
                      nodeicov("OT1") +
                      nodeicov("n1_indegree") +
                      edgecov(pt_secondarysim_norm) +
                      edgecov(pt_divisivesim_norm) +
                      edgecov(pt_interact_secon_divisive) +
                      edgecov(pt_n1) +
                      edges, 
                    eval.loglik = TRUE , check.degeneracy = TRUE , 
                    control = control.ergm ( seed = seed , MCMC.samplesize = 10000 , MCMC.interval = 1000, parallel = 11))

################################################ SWEDEN ##################################################################

# import data files

load("se_n1.RData")
load("se_n3.RData")
load("se_coresim_norm.RData")
load("se_secondarysim_norm.RData")
load("se_divisivesim_norm.RData")
load("se_interact_core_secon.RData")
load("se_interact_core_divisive.RData")
load("se_interact_secon_divisive.RData")

## exponential random graph models

SE_fit1 <- ergm(se_n3 ~  
                    gwesp(1, fixed = T) +
                    twopath +
                    gwodegree(1, fixed = T) +
                    mutual(same = NULL, diff = FALSE, by = NULL, keep = NULL) +
                    nodeicov("OT1") +
                    nodeicov("n1_indegree") +
                    edgecov(se_coresim_norm) +
                    edgecov(se_n1) +
                    edges, 
                  eval.loglik = TRUE , check.degeneracy = TRUE , 
                  control = control.ergm ( seed = seed , MCMC.samplesize = 10000 , MCMC.interval = 1000))

SE_fit2 <- ergm(se_n3 ~  
                     gwesp(1, fixed = T) +
                     twopath +
                     gwodegree(1, fixed = T) +
                     mutual(same = NULL, diff = FALSE, by = NULL, keep = NULL) +
                     nodeicov("OT1") +
                     nodeicov("n1_indegree") +
                     edgecov(se_secondarysim_norm) +
                     edgecov(se_n1) +
                     edges, 
                   eval.loglik = TRUE , check.degeneracy = TRUE , 
                   control = control.ergm ( seed = seed , MCMC.samplesize = 10000 , MCMC.interval = 1000))

SE_fit3 <- ergm(se_n3 ~
                     gwesp(1, fixed = T) +
                     twopath +
                     gwodegree(1, fixed = T) +
                     mutual(same = NULL, diff = FALSE, by = NULL, keep = NULL) +
                     nodeicov("OT1") +
                     nodeicov("n1_indegree") +
                     edgecov(se_divisivesim_norm) +
                     edgecov(se_n1) +
                     edges, 
                   eval.loglik = TRUE , check.degeneracy = TRUE , 
                   control = control.ergm ( seed = seed , MCMC.samplesize = 10000 , MCMC.interval = 1000))

## exponential random graph models with interactions

SE_fit4 <- ergm(se_n3 ~  
                      gwesp(1, fixed = T) +
                      twopath +
                      gwodegree(1, fixed = T) +
                      mutual(same = NULL, diff = FALSE, by = NULL, keep = NULL) +
                      nodeicov("OT1") +
                      nodeicov("n1_indegree") +
                      edgecov(se_coresim_norm) +
                      edgecov(se_secondarysim_norm) +
                      edgecov(se_interact_core_secon) +
                      edgecov(se_n1) +
                      edges, 
                    eval.loglik = TRUE , check.degeneracy = TRUE , 
                    control = control.ergm ( seed = seed , MCMC.samplesize = 10000 , MCMC.interval = 1000, parallel = 11))

SE_fit5 <- ergm(se_n3 ~  
                      gwesp(1, fixed = T) +
                      twopath +
                      gwodegree(1, fixed = T) +
                      mutual(same = NULL, diff = FALSE, by = NULL, keep = NULL) +
                      nodeicov("OT1") +
                      nodeicov("n1_indegree") +
                      edgecov(se_coresim_norm) +
                      edgecov(se_divisivesim_norm) +
                      edgecov(se_interact_core_divisive) +
                      edgecov(se_n1) +
                      edges, 
                    eval.loglik = TRUE , check.degeneracy = TRUE , 
                    control = control.ergm ( seed = seed , MCMC.samplesize = 10000 , MCMC.interval = 1000, parallel = 11))

SE_fit6 <- ergm(se_n3 ~  
                      gwesp(1, fixed = T) +
                      twopath +
                      gwodegree(1, fixed = T) +
                      mutual(same = NULL, diff = FALSE, by = NULL, keep = NULL) +
                      nodeicov("OT1") +
                      nodeicov("n1_indegree") +
                      edgecov(se_secondarysim_norm) +
                      edgecov(se_divisivesim_norm) +
                      edgecov(se_interact_secon_divisive) +
                      edgecov(se_n1) +
                      edges, 
                    eval.loglik = TRUE , check.degeneracy = TRUE , 
                    control = control.ergm ( seed = seed , MCMC.samplesize = 10000 , MCMC.interval = 1000, parallel = 11))

################################################## TAIWAN ################################################################

# import data files

load("tw_n1.RData")
load("tw_n3.RData")
load("tw_coresim_norm.RData")
load("tw_secondarysim_norm.RData")
load("tw_divisivesim_norm.RData")
load("tw_interact_core_secon.RData")
load("tw_interact_core_divisive.RData")
load("tw_interact_secon_divisive.RData")

## exponential random graph models

TW_fit1 <- ergm(tw_n3 ~ 
                    gwidegree(1, fixed = T) +
                    gwesp(1, fixed = T) +
                    twopath +
                    gwodegree(1, fixed = T) +
                    mutual(same = NULL, diff = FALSE, by = NULL, keep = NULL) +
                    nodeicov("OT1") +
                    nodeicov("n1_indegree") +
                    edgecov(tw_coresim_norm) +
                    edgecov(tw_n1) +
                    edges, 
                  eval.loglik = TRUE , check.degeneracy = TRUE , 
                  control = control.ergm ( seed = seed , MCMC.samplesize = 10000 , MCMC.interval = 1000))

TW_fit2 <- ergm(tw_n3 ~ 
                     gwidegree(1, fixed = T) +
                     gwesp(1, fixed = T) +
                     twopath +
                     gwodegree(1, fixed = T) +
                     mutual(same = NULL, diff = FALSE, by = NULL, keep = NULL) +
                     nodeicov("OT1") +
                     nodeicov("n1_indegree") +
                     edgecov(tw_secondarysim_norm) +
                     edgecov(tw_n1) +
                     edges, 
                   eval.loglik = TRUE , check.degeneracy = TRUE , 
                   control = control.ergm ( seed = seed , MCMC.samplesize = 10000 , MCMC.interval = 1000))

TW_fit3 <- ergm(tw_n3 ~ 
                     gwidegree(1, fixed = T) +
                     gwesp(1, fixed = T) +
                     twopath +
                     gwodegree(1, fixed = T) +
                     mutual(same = NULL, diff = FALSE, by = NULL, keep = NULL) +
                     nodeicov("OT1") +
                     nodeicov("n1_indegree") +
                     edgecov(tw_divisivesim_norm) +
                     edgecov(tw_n1) +
                     edges, 
                   eval.loglik = TRUE , check.degeneracy = TRUE , 
                   control = control.ergm ( seed = seed , MCMC.samplesize = 10000 , MCMC.interval = 1000))

## exponential random graph models with interactions

TW_fit4 <- ergm(tw_n3 ~ 
                      gwidegree(1, fixed = T) +
                      gwesp(1, fixed = T) +
                      twopath +
                      gwodegree(1, fixed = T) +
                      mutual(same = NULL, diff = FALSE, by = NULL, keep = NULL) +
                      nodeicov("OT1") +
                      nodeicov("n1_indegree") +
                      edgecov(tw_coresim_norm) +
                      edgecov(tw_secondarysim_norm) +
                      edgecov(tw_interact_core_secon) +
                      edgecov(tw_n1) +
                      edges, 
                    eval.loglik = TRUE , check.degeneracy = TRUE , 
                    control = control.ergm ( seed = seed , MCMC.samplesize = 10000 , MCMC.interval = 1000, parallel = 11))

TW_fit5 <- ergm(tw_n3 ~ 
                      gwidegree(1, fixed = T) +
                      gwesp(1, fixed = T) +
                      twopath +
                      gwodegree(1, fixed = T) +
                      mutual(same = NULL, diff = FALSE, by = NULL, keep = NULL) +
                      nodeicov("OT1") +
                      nodeicov("n1_indegree") +
                      edgecov(tw_coresim_norm) +
                      edgecov(tw_divisivesim_norm) +
                      edgecov(tw_interact_core_divisive) +
                      edgecov(tw_n1) +
                      edges, 
                    eval.loglik = TRUE , check.degeneracy = TRUE , 
                    control = control.ergm ( seed = seed , MCMC.samplesize = 10000 , MCMC.interval = 1000, parallel = 11))

TW_fit6 <- ergm(tw_n3 ~ 
                      gwidegree(1, fixed = T) +
                      gwesp(1, fixed = T) +
                      twopath +
                      gwodegree(1, fixed = T) +
                      mutual(same = NULL, diff = FALSE, by = NULL, keep = NULL) +
                      nodeicov("OT1") +
                      nodeicov("n1_indegree") +
                      edgecov(tw_secondarysim_norm) +
                      edgecov(tw_divisivesim_norm) +
                      edgecov(tw_interact_secon_divisive) +
                      edgecov(tw_n1) +
                      edges, 
                    eval.loglik = TRUE , check.degeneracy = TRUE , 
                    control = control.ergm ( seed = seed , MCMC.samplesize = 10000 , MCMC.interval = 1000, parallel = 11))


############################################## ERGM RESULT TABLES #########################################################

## ERGM RESULTS FOR MODELS WITHOUT INTERACTION

core_resuts_ergm <- screenreg(list(AU_fit1, BR_fit1, CZ_fit1, DE_fit1, FI_fit1, IE_fit1, JP_fit1, KR_fit1, PT_fit1, SE_fit1, TW_fit1), single.row = F, stars = c(0.001, 0.01, 0.05, 0.1))

secondary_resuts_ergm <- screenreg(list(AU_fit2, BR_fit2, CZ_fit2, DE_fit2, FI_fit2, IE_fit2, JP_fit2, KR_fit2, PT_fit2, SE_fit2, TW_fit2), single.row = F, stars = c(0.001, 0.01, 0.05, 0.1))

divisive_resuts_ergm <- screenreg(list(AU_fit3, BR_fit3, CZ_fit3, DE_fit3, FI_fit3, IE_fit3, JP_fit3, KR_fit3, PT_fit3, SE_fit3, TW_fit3), single.row = F, stars = c(0.001, 0.01, 0.05, 0.1))

## ERGM RESULTS FOR MODELS INCLUDIN INTERACTION TERMS

core_secon_interact_resuts_ergm <- screenreg(list(AU_fit4, BR_fit4, CZ_fit4, DE_fit4, FI_fit4, IE_fit4, JP_fit4, KR_fit4, PT_fit4, SE_fit4, TW_fit4), single.row = F, stars = c(0.001, 0.01, 0.05, 0.1))

core_divisive_interact_resuts_ergm <- screenreg(list(AU_fit5, BR_fit5, CZ_fit5, DE_fit5, FI_fit5, IE_fit5, JP_fit5, KR_fit5, PT_fit5, SE_fit5, TW_fit5), single.row = F, stars = c(0.001, 0.01, 0.05, 0.1))

secon_divisive_interact_resuts_ergm <- screenreg(list(AU_fit6, BR_fit6, CZ_fit6, DE_fit6, FI_fit6, IE_fit6, JP_fit6, KR_fit6, PT_fit6, SE_fit6, TW_fit6), single.row = F, stars = c(0.001, 0.01, 0.05, 0.1))

#############################################################################################################################

############################################# GOODNESS OF FIT ###############################################################

## gof plots for the best fitting model by country

AU_gofM1 <- gof(AU_fit1)

pdf("AU_gofM1.pdf")
par(mfrow=c(2,3))
par(oma=c(0.5,2,1,0.5))
plot(AU_gofM1, main = "Goodness of fit diagnostics: Australia M1")
dev.off()

BR_gofM3 <- gof(BR_fit3)

pdf("BR_gofM3.pdf")
par(mfrow=c(2,3))
par(oma=c(0.5,2,1,0.5))
plot(BR_gofM3, main = "Goodness of fit diagnostics: Brazil M3")
dev.off()

CZ_gofM3 <- gof(CZ_fit3)

pdf("CZ_gofM3.pdf")
par(mfrow=c(2,3))
par(oma=c(0.5,2,1,0.5))
plot(CZ_gofM3, main = "Goodness of fit diagnostics: the Czech Republic M3")
dev.off()

DE_gofM2 <- gof(DE_fit2)

pdf("DE_gofM2.pdf")
par(mfrow=c(2,3))
par(oma=c(0.5,2,1,0.5))
plot(DE_gofM2, main = "Goodness of fit diagnostics: Germany M2")
dev.off()

FI_gofM3 <- gof(FI_fit3)

pdf("FI_gofM3.pdf")
par(mfrow=c(2,3))
par(oma=c(0.5,2,1,0.5))
plot(FI_gofM3, main = "Goodness of fit diagnostics: Finland M3")
dev.off()

IE_gofM2 <- gof(IE_fit2)

pdf("IE_gofM2.pdf")
par(mfrow=c(2,3))
par(oma=c(0.5,2,1,0.5))
plot(IE_gofM2, main = "Goodness of fit diagnostics: Ireland M2")
dev.off()

JP_gofM3 <- gof(JP_fit3)

pdf("JP_gofM3.pdf")
par(mfrow=c(2,3))
par(oma=c(0.5,2,1,0.5))
plot(JP_gofM3, main = "Goodness of fit diagnostics: Japan M3")
dev.off()

KR_gofM3 <- gof(KR_fit3)

pdf("KR_gofM3.pdf")
par(mfrow=c(2,3))
par(oma=c(0.5,2,1,0.5))
plot(KR_gofM3, main = "Goodness of fit diagnostics: Korea M3")
dev.off()

PT_gofM2 <- gof(PT_fit2)

pdf("PT_gofM2.pdf")
par(mfrow=c(2,3))
par(oma=c(0.5,2,1,0.5))
plot(PT_gofM2, main = "Goodness of fit diagnostics: Portugal M2")
dev.off()

SE_gofM3 <- gof(SE_fit3)

pdf("SE_gofM3.pdf")
par(mfrow=c(2,3))
par(oma=c(0.5,2,1,0.5))
plot(SE_gofM3, main = "Goodness of fit diagnostics: Sweden M3")
dev.off()

TW_gofM1 <- gof(TW_fit1)

pdf("TW_gofM1.pdf")
par(mfrow=c(2,3))
par(oma=c(0.5,2,1,0.5))
plot(TW_gofM1, main = "Goodness of fit diagnostics: Taiwan M1")
dev.off()

#############################################################################################################################

## MCMC diagnostics for the best fitting model by country

## AUSTRALIA

# Loading object
object_AU <- AU_fit1
# Extract MCMC sample
sm_AU <- NVL3(object_AU[["sample"]], as.mcmc.list(.))

# plot MCMC
pdf("AU_mcmc_M1.pdf", width = 21 , height = 30)
par(mfrow=c(5,2),oma = c(0, 0, 2, 0))
for(i in 1:ncol(sm_AU[[1]])){
  # Setting values
  yvals <- sm_AU[[1]][,i]
  index.vals <- 1:length(sm_AU[[1]][,i])
  # Plotting
  plot(index.vals, yvals, type = 'l', main = paste(colnames(sm_AU[[1]])[i]))
  plot(density(yvals), main = paste(colnames(sm_AU[[1]])[i]))
}
mtext("MCMC diagnostics: Australia M1", outer = TRUE, cex = 1.5)
dev.off()

## BRAZIL

# Loading object
object_BR <- BR_fit3
# Extract MCMC sample
sm_BR <- NVL3(object_BR[["sample"]], as.mcmc.list(.))

# plot MCMC
pdf("BR_mcmc_M3.pdf", width = 21 , height = 30)
par(mfrow=c(5,2),oma = c(0, 0, 2, 0))
for(i in 1:ncol(sm_BR[[1]])){
  # Setting values
  yvals <- sm_BR[[1]][,i]
  index.vals <- 1:length(sm_BR[[1]][,i])
  # Plotting
  plot(index.vals, yvals, type = 'l', main = paste(colnames(sm_BR[[1]])[i]))
  plot(density(yvals), main = paste(colnames(sm_BR[[1]])[i]))
}
mtext("MCMC diagnostics: Brazil M3", outer = TRUE, cex = 1.5)
dev.off()

## THE CZECH REPUBLIC

# Loading object
object_CZ <- CZ_fit3
# Extract MCMC sample
sm_CZ <- NVL3(object_CZ[["sample"]], as.mcmc.list(.))

# plot MCMC
pdf("CZ_mcmc_M3.pdf", width = 21 , height = 30)
par(mfrow=c(5,2),oma = c(0, 0, 2, 0))
for(i in 1:ncol(sm_CZ[[1]])){
  # Setting values
  yvals <- sm_CZ[[1]][,i]
  index.vals <- 1:length(sm_CZ[[1]][,i])
  # Plotting
  plot(index.vals, yvals, type = 'l', main = paste(colnames(sm_CZ[[1]])[i]))
  plot(density(yvals), main = paste(colnames(sm_CZ[[1]])[i]))
}
mtext("MCMC diagnostics: the Czech Republic M3", outer = TRUE, cex = 1.5)
dev.off()

## GERMANY

# Loading object
object_DE <- DE_fit2
# Extract MCMC sample
sm_DE <- NVL3(object_DE[["sample"]], as.mcmc.list(.))

# plot MCMC
pdf("DE_mcmc_M2.pdf", width = 21 , height = 30)
par(mfrow=c(5,2),oma = c(0, 0, 2, 0))
for(i in 1:ncol(sm_DE[[1]])){
  # Setting values
  yvals <- sm_DE[[1]][,i]
  index.vals <- 1:length(sm_DE[[1]][,i])
  # Plotting
  plot(index.vals, yvals, type = 'l', main = paste(colnames(sm_DE[[1]])[i]))
  plot(density(yvals), main = paste(colnames(sm_DE[[1]])[i]))
}
mtext("MCMC diagnostics: Germany M2", outer = TRUE, cex = 1.5)
dev.off()

## FINLAND

# Loading object
object_FI <- FI_fit3
# Extract MCMC sample
sm_FI <- NVL3(object_FI[["sample"]], as.mcmc.list(.))

# plot MCMC
pdf("FI_mcmc_M3.pdf", width = 21 , height = 30)
par(mfrow=c(5,2),oma = c(0, 0, 2, 0))
for(i in 1:ncol(sm_FI[[1]])){
  # Setting values
  yvals <- sm_FI[[1]][,i]
  index.vals <- 1:length(sm_FI[[1]][,i])
  # Plotting
  plot(index.vals, yvals, type = 'l', main = paste(colnames(sm_FI[[1]])[i]))
  plot(density(yvals), main = paste(colnames(sm_FI[[1]])[i]))
}
mtext("MCMC diagnostics: Finland M3", outer = TRUE, cex = 1.5)
dev.off()

## IRELAND

# Loading object
object_IE <- IE_fit2
# Extract MCMC sample
sm_IE <- NVL3(object_IE[["sample"]], as.mcmc.list(.))

# plot MCMC
pdf("IE_mcmc_M2.pdf", width = 21 , height = 30)
par(mfrow=c(5,2),oma = c(0, 0, 2, 0))
for(i in 1:ncol(sm_IE[[1]])){
  # Setting values
  yvals <- sm_IE[[1]][,i]
  index.vals <- 1:length(sm_IE[[1]][,i])
  # Plotting
  plot(index.vals, yvals, type = 'l', main = paste(colnames(sm_IE[[1]])[i]))
  plot(density(yvals), main = paste(colnames(sm_IE[[1]])[i]))
}
mtext("MCMC diagnostics: Ireland M2", outer = TRUE, cex = 1.5)
dev.off()

## JAPAN

# Loading object
object_JP <- JP_fit3
# Extract MCMC sample
sm_JP <- NVL3(object_JP[["sample"]], as.mcmc.list(.))

# plot MCMC
pdf("JP_mcmc_M3.pdf", width = 21 , height = 30)
par(mfrow=c(5,2),oma = c(0, 0, 2, 0))
for(i in 1:ncol(sm_JP[[1]])){
  # Setting values
  yvals <- sm_JP[[1]][,i]
  index.vals <- 1:length(sm_JP[[1]][,i])
  # Plotting
  plot(index.vals, yvals, type = 'l', main = paste(colnames(sm_JP[[1]])[i]))
  plot(density(yvals), main = paste(colnames(sm_JP[[1]])[i]))
}
mtext("MCMC diagnostics: Japan M3", outer = TRUE, cex = 1.5)
dev.off()

## KOREA

# Loading object
object_KR <- KR_fit3
# Extract MCMC sample
sm_KR <- NVL3(object_KR[["sample"]], as.mcmc.list(.))

# plot MCMC
pdf("KR_mcmc_M3.pdf", width = 21 , height = 30)
par(mfrow=c(5,2),oma = c(0, 0, 2, 0))
for(i in 1:ncol(sm_KR[[1]])){
  # Setting values
  yvals <- sm_KR[[1]][,i]
  index.vals <- 1:length(sm_KR[[1]][,i])
  # Plotting
  plot(index.vals, yvals, type = 'l', main = paste(colnames(sm_KR[[1]])[i]))
  plot(density(yvals), main = paste(colnames(sm_KR[[1]])[i]))
}
mtext("MCMC diagnostics: Korea M3", outer = TRUE, cex = 1.5)
dev.off()

## PORTUGAL

# Loading object
object_PT <- PT_fit2
# Extract MCMC sample
sm_PT <- NVL3(object_PT[["sample"]], as.mcmc.list(.))

# plot MCMC
pdf("PT_mcmc_M2.pdf", width = 21 , height = 30)
par(mfrow=c(5,2),oma = c(0, 0, 2, 0))
for(i in 1:ncol(sm_PT[[1]])){
  # Setting values
  yvals <- sm_PT[[1]][,i]
  index.vals <- 1:length(sm_PT[[1]][,i])
  # Plotting
  plot(index.vals, yvals, type = 'l', main = paste(colnames(sm_PT[[1]])[i]))
  plot(density(yvals), main = paste(colnames(sm_PT[[1]])[i]))
}
mtext("MCMC diagnostics: Portugal M2", outer = TRUE, cex = 1.5)
dev.off()

## SWEDEN

# Loading object
object_SE <- SE_fit3
# Extract MCMC sample
sm_SE <- NVL3(object_SE[["sample"]], as.mcmc.list(.))

# plot MCMC
pdf("SE_mcmc_M3.pdf", width = 21 , height = 30)
par(mfrow=c(5,2),oma = c(0, 0, 2, 0))
for(i in 1:ncol(sm_SE[[1]])){
  # Setting values
  yvals <- sm_SE[[1]][,i]
  index.vals <- 1:length(sm_SE[[1]][,i])
  # Plotting
  plot(index.vals, yvals, type = 'l', main = paste(colnames(sm_SE[[1]])[i]))
  plot(density(yvals), main = paste(colnames(sm_SE[[1]])[i]))
}
mtext("MCMC diagnostics: Sweden M3", outer = TRUE, cex = 1.5)
dev.off()

## TAIWAN

# Loading object
object_TW <- TW_fit1
# Extract MCMC sample
sm_TW <- NVL3(object_TW[["sample"]], as.mcmc.list(.))

# plot MCMC
pdf("TW_mcmc_M1.pdf", width = 21 , height = 30)
par(mfrow=c(5,2),oma = c(0, 0, 2, 0))
for(i in 1:ncol(sm_TW[[1]])){
  # Setting values
  yvals <- sm_TW[[1]][,i]
  index.vals <- 1:length(sm_TW[[1]][,i])
  # Plotting
  plot(index.vals, yvals, type = 'l', main = paste(colnames(sm_TW[[1]])[i]))
  plot(density(yvals), main = paste(colnames(sm_TW[[1]])[i]))
}
mtext("MCMC diagnostics: Taiwan M1", outer = TRUE, cex = 1.5)
dev.off()

#################################################################

## MARGINAL EFFECT PLOTS AND EDGE PROBABILITIES (PLOTS)
## CODE APPLIED FROM:
# Heaney, Michael T., and Philip Leifeld. 2018. 
# "Contributions by Interest Groups to Lobbying Coalitions." 
# The Journal of Politics 80 (2): 494–509.
# doi: 10.1086/694545. 
# https://www.journals.uchicago.edu/doi/abs/10.1086/694545.

#################################################################

## required packages

if(!require(network)) install.packages("inline",repos = "http://cran.us.r-project.org")
library("network") 
if(!require(sna)) install.packages("inline",repos = "http://cran.us.r-project.org")
library("sna") 
if(!require(ergm)) install.packages("inline",repos = "http://cran.us.r-project.org")
library("ergm") 
if(!require(inline)) install.packages("inline",repos = "http://cran.us.r-project.org")
library("inline") 
if(!require(Rcpp)) install.packages("inline",repos = "http://cran.us.r-project.org")
library("Rcpp") 
if(!require(btergm)) install.packages("btergm",repos = "http://cran.us.r-project.org")
library("btergm")
if(!require(ggplot2)) install.packages("ggplot2",repos = "http://cran.us.r-project.org")
library("ggplot2") 
if(!require(reshape2)) install.packages("reshape2",repos = "https://github.com/hadley/reshape")
library("reshape2") 
if(!require(grid)) install.packages("inline",repos = "http://cran.us.r-project.org")
library("grid")
if(!require(gridExtra)) install.packages("inline",repos = "http://cran.us.r-project.org")
library("gridExtra")

################################# PREDICTED PROPABILITIES ############################################

## function for plotting facets with variable 1 conditional on variable 2 
facets <- function(edgeprobs, country, var1, var2, varname1, varname2) {   
  
  dyads <- edgeprobs  
  
  # cut network embeddedness into slices   
  dyads$v1 <- c(dyads[var1])[[1]]   
  v2 <- c(dyads[var2])[[1]]   
  v2.quantiles <- quantile(v2)   
  v2 <- cut(v2, v2.quantiles, labels = names(v2.quantiles)[2:5])   
  v2[is.na(v2)] <- "25%"   
  dta <- transform(dyads, v2 = v2)   
  
  # plot conditional probabilities with facets   
  gp <- ggplot(data = dta, aes(x = v1, y = probability))   
  print(gp + stat_smooth(method = "lm", fullrange = TRUE, color = "black") + 
          facet_wrap( ~ v2,       
                      ncol = 2) + xlab(varname1) + ylab("Probability") + ggtitle(
                        paste0(country, ": effect of ", varname1, " conditional on ", 
                               varname2)))   
} 

############## AUSTRALIA ############################

# create dyadic datasets 
edgeprob.au_cdi <- edgeprob(AU_fit5) 

facets(edgeprobs = edgeprob.au_cdi, country = "Australia",         
       var1 = "edgecov.au_divisivesim_norm[[i]]",         
       var2 = "edgecov.au_coresim_norm[[i]]",         
       varname1 = "divisive belief similarity", 
       varname2 = "core belief similarity") 

# marginal effects plots
marginalplot(AU_fit5, 
             var1 = "edgecov.au_divisivesim_norm", 
             var2 = "edgecov.au_coresim_norm",
             inter = "edgecov.au_interact_core_divisive", 
             ylab = "divisive belief similarity",               
             xlab = "core belief similarity",               
             rug = TRUE) + 
  ggtitle("AU core-divisive interaction") + theme_bw()

############## BRAZIL ############################

# create dyadic datasets 
edgeprob.br_cdi <- edgeprob(BR_fit5) 

facets(edgeprobs = edgeprob.br_cdi, country = "Brazil",         
       var1 = "edgecov.br_divisivesim_norm[[i]]",         
       var2 = "edgecov.br_coresim_norm[[i]]",         
       varname1 = "divisive belief similarity", 
       varname2 = "core belief similarity") 

# marginal effects plots
marginalplot(BR_fit5, 
             var1 = "edgecov.br_divisivesim_norm", 
             var2 = "edgecov.br_coresim_norm",
             inter = "edgecov.br_interact_core_divisive", 
             ylab = "divisive belief similarity",               
             xlab = "core belief similarity",               
             rug = TRUE) + 
  ggtitle("BR core-divisive interaction") + theme_bw()

############## the CZECH REPUBLIC ############################

# create dyadic datasets 
edgeprob.cz_cdi <- edgeprob(CZ_fit5) 

facets(edgeprobs = edgeprob.cz_cdi, country = "the Czech Republic",         
       var1 = "edgecov.cz_divisivesim_norm[[i]]",         
       var2 = "edgecov.cz_coresim_norm[[i]]",         
       varname1 = "divisive belief similarity", 
       varname2 = "core belief similarity") 

# marginal effects plots
marginalplot(CZ_fit5, 
             var1 = "edgecov.cz_divisivesim_norm", 
             var2 = "edgecov.cz_coresim_norm",
             inter = "edgecov.cz_interact_core_divisive", 
             ylab = "divisive belief similarity",               
             xlab = "core belief similarity",               
             rug = TRUE) + 
  ggtitle("CZ core-divisive interaction") + theme_bw()

############## GERMANY ############################

# create dyadic datasets 
edgeprob.de_cdi <- edgeprob(DE_fit5) 

facets(edgeprobs = edgeprob.de_cdi, country = "Germany",         
       var1 = "edgecov.de_divisivesim_norm[[i]]",         
       var2 = "edgecov.de_coresim_norm[[i]]",         
       varname1 = "divisive belief similarity", 
       varname2 = "core belief similarity") 

# marginal effects plots
marginalplot(DE_fit5, 
             var1 = "edgecov.de_divisivesim_norm", 
             var2 = "edgecov.de_coresim_norm",
             inter = "edgecov.de_interact_core_divisive", 
             ylab = "divisive belief similarity",               
             xlab = "core belief similarity",               
             rug = TRUE) + 
  ggtitle("DE core-divisive interaction") + theme_bw()

############## FINLAND ############################

# create dyadic datasets 
edgeprob.fi_cdi <- edgeprob(FI_fit5) 

facets(edgeprobs = edgeprob.fi_cdi, country = "Finland",         
       var1 = "edgecov.fi_divisivesim_norm[[i]]",         
       var2 = "edgecov.fi_coresim_norm[[i]]",         
       varname1 = "divisive belief similarity", 
       varname2 = "core belief similarity") 


# marginal effects plots
marginalplot(FI_fit5, 
             var1 = "edgecov.fi_divisivesim_norm", 
             var2 = "edgecov.fi_coresim_norm",
             inter = "edgecov.fi_interact_core_divisive", 
             ylab = "divisive belief similarity",               
             xlab = "core belief similarity",               
             rug = TRUE) + 
  ggtitle("FI core-divisive interaction") + theme_bw()

############## IRELAND ############################

# create dyadic datasets 
edgeprob.ie_cdi <- edgeprob(IE_fit5) 

facets(edgeprobs = edgeprob.ie_cdi, country = "Ireland",         
       var1 = "edgecov.ie_divisivesim_norm[[i]]",         
       var2 = "edgecov.ie_coresim_norm[[i]]",         
       varname1 = "divisive belief similarity", 
       varname2 = "core belief similarity") 

# marginal effects plots
marginalplot(IE_fit5, 
             var1 = "edgecov.ie_divisivesim_norm", 
             var2 = "edgecov.ie_coresim_norm",
             inter = "edgecov.ie_interact_core_divisive", 
             ylab = "divisive belief similarity",               
             xlab = "core belief similarity",               
             rug = TRUE) + 
  ggtitle("IE core-divisive interaction") + theme_bw()

############## JAPAN ############################

# create dyadic datasets 
edgeprob.jp_cdi <- edgeprob(JP_fit5) 

facets(edgeprobs = edgeprob.jp_cdi, country = "Japan",         
       var1 = "edgecov.jp_divisivesim_norm[[i]]",         
       var2 = "edgecov.jp_coresim_norm[[i]]",         
       varname1 = "divisive belief similarity", 
       varname2 = "core belief similarity") 

# marginal effects plots
marginalplot(JP_fit5, 
             var1 = "edgecov.jp_divisivesim_norm", 
             var2 = "edgecov.jp_coresim_norm",
             inter = "edgecov.jp_interact_core_divisive", 
             ylab = "divisive belief similarity",               
             xlab = "core belief similarity",               
             rug = TRUE) + 
  ggtitle("JP core-divisive interaction") + theme_bw()

############## KOREA ############################

# create dyadic datasets 
edgeprob.kr_cdi <- edgeprob(KR_fit5) 

facets(edgeprobs = edgeprob.kr_cdi, country = "Korea",         
       var1 = "edgecov.kr_divisivesim_norm[[i]]",         
       var2 = "edgecov.kr_coresim_norm[[i]]",         
       varname1 = "divisive belief similarity", 
       varname2 = "core belief similarity") 


# marginal effects plots
marginalplot(KR_fit5, 
             var1 = "edgecov.kr_divisivesim_norm", 
             var2 = "edgecov.kr_coresim_norm",
             inter = "edgecov.kr_interact_core_divisive", 
             ylab = "divisive belief similarity",               
             xlab = "core belief similarity",               
             rug = TRUE) + 
  ggtitle("KR core-divisive interaction") + theme_bw()

############## PORTUGAL ############################

# create dyadic datasets 
edgeprob.pt_cdi <- edgeprob(PT_fit5) 

facets(edgeprobs = edgeprob.pt_cdi, country = "Portugal",         
       var1 = "edgecov.pt_divisivesim_norm[[i]]",         
       var2 = "edgecov.pt_coresim_norm[[i]]",         
       varname1 = "divisive belief similarity", 
       varname2 = "core belief similarity") 


# marginal effects plots
marginalplot(PT_fit5, 
             var1 = "edgecov.pt_divisivesim_norm", 
             var2 = "edgecov.pt_coresim_norm",
             inter = "edgecov.pt_interact_core_divisive", 
             ylab = "divisive belief similarity",               
             xlab = "core belief similarity",               
             rug = TRUE) + 
  ggtitle("PT core-divisive interaction") + theme_bw()

############## SWEDEN ############################

# create dyadic datasets 
edgeprob.se_cdi <- edgeprob(SE_fit5) 

facets(edgeprobs = edgeprob.se_cdi, country = "Sweden",         
       var1 = "edgecov.se_divisivesim_norm[[i]]",         
       var2 = "edgecov.se_coresim_norm[[i]]",         
       varname1 = "divisive belief similarity", 
       varname2 = "core belief similarity") 

# marginal effects plots
marginalplot(SE_fit5, 
             var1 = "edgecov.se_divisivesim_norm", 
             var2 = "edgecov.se_coresim_norm",
             inter = "edgecov.se_interact_core_divisive", 
             ylab = "divisive belief similarity",               
             xlab = "core belief similarity",               
             rug = TRUE) + 
  ggtitle("SE core-divisive interaction") + theme_bw()

############## TAIWAN ############################

# create dyadic datasets 
edgeprob.tw_cdi <- edgeprob(TW_fit5) 

facets(edgeprobs = edgeprob.tw_cdi, country = "Taiwan",         
       var1 = "edgecov.tw_divisivesim_norm[[i]]",         
       var2 = "edgecov.tw_coresim_norm[[i]]",         
       varname1 = "divisive belief similarity", 
       varname2 = "core belief similarity") 


# marginal effects plots
marginalplot(TW_fit5, 
             var1 = "edgecov.tw_divisivesim_norm", 
             var2 = "edgecov.tw_coresim_norm",
             inter = "edgecov.tw_interact_core_divisive", 
             ylab = "divisive belief similarity",               
             xlab = "core belief similarity",               
             rug = TRUE) + 
  ggtitle("TW core-divisive interaction") + theme_bw()

################# MULTIPLE PLOTS/PAGE ############

# function for plotting facets with variable 1 conditional on variable 2, for plots with one title for multiple countries
facets2 <- function(edgeprobs, country, var1, var2, varname1, varname2) {   
  
  dyads <- edgeprobs  
  
  # cut network embeddedness into slices   
  dyads$v1 <- c(dyads[var1])[[1]]   
  v2 <- c(dyads[var2])[[1]]   
  v2.quantiles <- quantile(v2)   
  v2 <- cut(v2, v2.quantiles, labels = names(v2.quantiles)[2:5])   
  v2[is.na(v2)] <- "25%"   
  dta <- transform(dyads, v2 = v2)   
  
  # plot conditional probabilities with facets   
  gp <- ggplot(data = dta, aes(x = v1, y = probability))   
  print(gp + stat_smooth(method = "lm", fullrange = TRUE, color = "black") + 
          facet_wrap( ~ v2,       
                      ncol = 2) + xlab(varname1) + ylab("Probability") + ggtitle(
                        paste0(country)))   
} 

e1 <- facets2(edgeprobs = edgeprob.au_cdi, country = "Australia",         
              var1 = "edgecov.au_divisivesim_norm[[i]]",         
              var2 = "edgecov.au_coresim_norm[[i]]",         
              varname1 = "divisive belief similarity", 
              varname2 = "core belief similarity") 

e2 <- facets2(edgeprobs = edgeprob.br_cdi, country = "Brazil",         
              var1 = "edgecov.br_divisivesim_norm[[i]]",         
              var2 = "edgecov.br_coresim_norm[[i]]",         
              varname1 = "divisive belief similarity", 
              varname2 = "core belief similarity") 

e3 <- facets2(edgeprobs = edgeprob.cz_cdi, country = "the Czech Republic",         
              var1 = "edgecov.cz_divisivesim_norm[[i]]",         
              var2 = "edgecov.cz_coresim_norm[[i]]",         
              varname1 = "divisive belief similarity", 
              varname2 = "core belief similarity") 

e4 <- facets2(edgeprobs = edgeprob.de_cdi, country = "Germany",         
              var1 = "edgecov.de_divisivesim_norm[[i]]",         
              var2 = "edgecov.de_coresim_norm[[i]]",         
              varname1 = "divisive belief similarity", 
              varname2 = "core belief similarity") 

e5 <- facets2(edgeprobs = edgeprob.fi_cdi, country = "Finland",         
              var1 = "edgecov.fi_divisivesim_norm[[i]]",         
              var2 = "edgecov.fi_coresim_norm[[i]]",         
              varname1 = "divisive belief similarity", 
              varname2 = "core belief similarity") 

e6 <- facets2(edgeprobs = edgeprob.ie_cdi, country = "Ireland",         
              var1 = "edgecov.ie_divisivesim_norm[[i]]",         
              var2 = "edgecov.ie_coresim_norm[[i]]",         
              varname1 = "divisive belief similarity", 
              varname2 = "core belief similarity") 

e7 <- facets2(edgeprobs = edgeprob.jp_cdi, country = "Japan",         
              var1 = "edgecov.jp_divisivesim_norm[[i]]",         
              var2 = "edgecov.jp_coresim_norm[[i]]",         
              varname1 = "divisive belief similarity", 
              varname2 = "core belief similarity") 

e8 <- facets2(edgeprobs = edgeprob.kr_cdi, country = "Korea",         
              var1 = "edgecov.kr_divisivesim_norm[[i]]",         
              var2 = "edgecov.kr_coresim_norm[[i]]",         
              varname1 = "divisive belief similarity", 
              varname2 = "core belief similarity") 

e9 <- facets2(edgeprobs = edgeprob.pt_cdi, country = "Portugal",         
              var1 = "edgecov.pt_divisivesim_norm[[i]]",         
              var2 = "edgecov.pt_coresim_norm[[i]]",         
              varname1 = "divisive belief similarity", 
              varname2 = "core belief similarity") 

e10 <- facets2(edgeprobs = edgeprob.se_cdi, country = "Sweden",         
               var1 = "edgecov.se_divisivesim_norm[[i]]",         
               var2 = "edgecov.se_coresim_norm[[i]]",         
               varname1 = "divisive belief similarity", 
               varname2 = "core belief similarity") 

e11 <- facets2(edgeprobs = edgeprob.tw_cdi, country = "Taiwan",         
               var1 = "edgecov.tw_divisivesim_norm[[i]]",         
               var2 = "edgecov.tw_coresim_norm[[i]]",         
               varname1 = "divisive belief similarity", 
               varname2 = "core belief similarity") 


pdf("edgeprobs_appendix1.pdf")
grid.newpage()
pushViewport(viewport(layout = grid.layout(3, 2, heights = unit(c(0.5, 5, 5), "null"))))   
vplayout <- function(x, y) viewport(layout.pos.row = x, layout.pos.col = y)
print(e1, vp = vplayout(2, 1))
print(e2, vp = vplayout(2, 2))
print(e3, vp = vplayout(3, 1))
print(e4, vp = vplayout(3, 2))
grid.text("Effect of divisive belief similarity on collaboration, conditional on policy core belief similarity", vp = viewport(layout.pos.row = 1, layout.pos.col = 1:2))
dev.off()

pdf("edgeprobs_appendix2.pdf")
grid.newpage()
pushViewport(viewport(layout = grid.layout(3, 2, heights = unit(c(0.5, 5, 5), "null"))))   
vplayout <- function(x, y) viewport(layout.pos.row = x, layout.pos.col = y)
print(e5, vp = vplayout(2, 1))
print(e6, vp = vplayout(2, 2))
print(e7, vp = vplayout(3, 1))
print(e8, vp = vplayout(3, 2))
grid.text("Effect of divisive belief similarity on collaboration, conditional on core belief similarity", vp = viewport(layout.pos.row = 1, layout.pos.col = 1:2))
dev.off()

pdf("edgeprobs_appendix3.pdf")
grid.newpage()
pushViewport(viewport(layout = grid.layout(3, 2, heights = unit(c(0.5, 5, 5), "null"))))   
vplayout <- function(x, y) viewport(layout.pos.row = x, layout.pos.col = y)
print(e9, vp = vplayout(2, 1))
print(e10, vp = vplayout(2, 2))
print(e11, vp = vplayout(3, 1))
grid.text("Effect of divisive belief similarity on collaboration, conditional on core belief similarity", vp = viewport(layout.pos.row = 1, layout.pos.col = 1:2))
dev.off()

